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Abstract. The standard pair approximation equations (PA) for the 
Susceptible-Infective-Recovered-Susceptible (SIRS) model of infection 
spread on a network of homogeneous degree k predict a thin phase of 
sustained oscillations for parameter values that correspond to diseases 
that confer long lasting immunity. Here we present a study of the depen- 
dence of this oscillatory phase on the parameter k and of its relevance 
to understand the behaviour of simulations on networks. For k = 4, we 
compare the phase diagram of the PA model with the results of simula- 
tions on regular random graphs (RRG) of the same degree. We show that 
for parameter values in the oscillatory phase, and even for large system 
sizes, the simulations either die out or exhibit damped oscillations, de- 
pending on the initial conditions. This failure of the standard PA model 
to capture the qualitative behaviour of the simulations on large RRGs is 
currently being investigated. 
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A number of approaches has been used to study the spreading dynamics of an 
infectious disease. A common paradigm, emerging from a simple deterministic 
framework, is to assume that populations are not spatially distributed so that 
individuals mix perfectly and contact each other with equal probability. Thus in 
the limit of infinite populations, the time evolution of the disease is described 
in terms of the densities of infectives and susceptibles as a function of time, and 
governed by a system of coupled ordinary differential equations which can be 
deduced from the law of mass action [Tj . Another approach is to use stochastic 
dynamics on a lattice (or more general graphs) where the variables at each node 
represent the state of an individual. The effects of spatial correlations that mass 
action models disregard play an important role in the behaviour of infection 
dynamics on graphs, and therefore also in real populations [2J. The ordinary 
pair approximation (PA) as well as various improvements to include higher order 
correlations have been proposed in the context of ecological and epidemiological 
deterministic models [3]. In 0], the performance of the PA in the description 
of the steady states and the dynamics of the Susceptible-Infective-Recovered- 
Susceptible (SIRS) model on the hypercubic lattice was analyzed in detail. 
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In this study, we consider the dynamics of the same epidemic model on a 
random network of homogeneous degree k and N nodes, a regular random graph 
of degree k (RRG-/c) . Each node can be occupied by an individual in susceptible 
(S), infected (I), or recovered (R) state. Infected individuals recover at rate S, 
recovered individuals lose immunity at rate 7, and infection of the susceptible 
node occurs at infection rate A multiplied by the number of its infected nearest 
neighbours n, n £ {0, 1, . . . , k}: 

I±R, 

R^S, (1) 
S^I . 

In the infinite population limit, with the assumptions of spatial homogeneity 
and uncorrelated pairs, the system is described by the deterministic equations 
of the standard or uncorrelated PA [3]: 
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In the above equations the variables s, i stand for the probability that a randomly 
chosen node is in state 5, /, and the variables si, sr, ri stand for the probability 
that a randomly chosen pair of nearest neighbour nodes is an SI, SR, RI pair. As 
expected, neglecting the pair correlations and setting the pair state probabilities 
equal to the product of the node state probabilities these equations reduce to 
the classic equations of the randomly mixed SIRS model. 

The phase diagram of the PA SIRS model ^ for k = 4 is plotted in Fig. 
QJ,). We have set the time scale so that 5 = 1. Region I represents susceptible- 
absorbing states and region II corresponds to active states that can be asymp- 
totically stable nodes or asymptotically stable foci. The critical line separating 
regions I and II corresponds to the transcritical bifurcation curve that is given 
by A c (7) = (7 + l)/(37 + 2) (black dashed line). In addition, for small values 
of 7 we find a new phase boundary (black solid line), that corresponds to a 
supercritical Hopf bifurcation of the nontrivial equilibrium and has been missed 
in previous studies of this model [4] . This boundary separates the active phase 
with constant densities from an active phase with oscillatory behavior, that is 
stable at low 7. 

In the thin phase of region III, the PA model predicts sustained oscillations 
in the thermodynamic limit. We have performed a systematic study of the de- 
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Fig. 1. a) Phase diagram in the (A, 7) plane for the PA SIRS model, where region I 
corresponds to susceptible-absorbing states and region II corresponds to active states 
with nonzero infective densities. The critical line between regions I and II is the black 
dashed line. The second critical curve (black solid line) bounds a region with limit cycle 
solutions (region III). Parameters: 5 = 1, k = 4. b) Phase diagram of the PA SIRS 
model for 8 — 1 and k — 2.1 (blue), k = 3 (red), k = 4 (black), k — 5 (green). Dashed 
(dotted) lines correspond to the transcritical (supercritical Hopf) bifurcation curves. 




pendence of this oscillatory phase on the parameter k and of its relevance to 
understand the behaviour of simulations on networks. 

The phase diagram of the PA model for 6=1 and several values of k in the 
range k > 2 is shown in Fig. [TJa). The critical lines separating the absorbing 
and the active phases (dashed lines) are given by A c (7) = (7 + l)/((fc — 1)7 + 
k — 2). Within the active phase, the dotted lines are numerical plots of Hopf 
bifurcation curves. The oscillatory phase is large for k > 2 and it gets thinner 
as k increases, but it persists for the whole range of 2 < k < 6. A similar phase 
diagram, with a Hopf bifurcation critical line bounding an oscillatory phase, 
was reported in other studies of related models [5], where SIR dynamics with 
different mechanisms of replenishment of susceptibles is modelled at the level 
of pairs with the standard or another closure approximation. These different 
models all exhibit an oscillatory phase in the regime of slow driving through 
introduction of new susceptible individuals (small 7 in the present case). This 
suggests that this oscillatory phase may be related with the phenomenon of 
recurrent epidemics in infectious diseases that confer permanent or long lasting 
immunity. 

We have compared the behaviour of the PA SIRS model © for k = 4 with 
the results of stochastic simulations on RRG-4 for several system sizes. In the 
stochastic simulations, the system was set in a random initial condition with 
given node and pair densities and an efficient algorithm for stochastic processes in 
spatially structured systems based on Gillespie's method [6] was used to update 
the states of the nodes according to the processes of infection, recovery and 
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immunity waning JTJ). For each set of parameter values and initial conditions, 
the simulations were averaged over 10 3 realizations of the RRG-4 graph. 
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Fig. 2. For 5 = 1, k = 4, comparison of the solutions of the PA deterministic model 
(dashed lines) with the results of stochastic simulations (solid lines) on a RRG-4 with 
N — 10 6 for parameter values in region II. Susceptible (infective) densities are plotted 
in blue (black). Parameters: a) 7 = 2.5, A = 2.5; b) 7 = 0.1, A = 2.5. 

The results of stochastic simulations for N = 10 6 and solutions of the PA 
SIRS equations ^ are shown in Figs. [2] and [3] The susceptible (blue lines) and 
the infective (black lines) densities are shown in Fig. [5] for two sets of parameter 
values: 7 = 2.5, A = 2.5 (Fig.[2t)) and 7 = 0.1, A = 2.5 (Fig.Gb)). The numerical 
solutions of the PA SIRS equations are plotted in dashed lines, and the results 
of the simulations in solid lines. For parameter values well within region II of the 
phase diagram as in Fig. [5^) there is excellent agreement between the solutions of 
the PA SIRS model for the same initial densities and the results of the stochastic 
simulations, both for the transient behaviour and for the steady states. This 
agreement deteriorates as 7 decreases and the boundary of the oscillatory region 
is approached as can be seen in Fig. [2b) ■ For parameter values in the oscillatory 
region III most simulations (black solid line) die out after a short transient (Fig. 
Eh>)) while the corresponding solutions of the PA SIRS deterministic model (blue 
solid line) converge to the stable limit cycle for all initial conditions (a typical set 
is denoted by B in the plot). By choosing initial conditions not far from the stable 
cycle predicted by the PA SIRS model to avoid extreme susceptible depletion 
during the transient, damped oscillations towards a non trivial equilibrium may 
also be observed in region III. In Fig. [3^,) a plot is shown of one of these surviving 
simulations (black solid line), together with the solution of the PA equations 
(blue solid line) for the same parameter values and initial conditions (in the plot 
denoted by A). Thus, instead of an oscillatory phase, the stochastic model on 
RRGs exhibits in region III a bistability phase, even for large system sizes. 

This failure of the PA model to capture the qualitative behaviour of the 
simulations on large RRGs is currently being investigated. Extinctions due to 
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finite size are one of the reasons why the oscillatory phase is seen as an ab- 
sorbing phase in the stochastic simulations. Indeed, as can be seen in Fig. [3J 
the oscillations predicted by the PA SIRS deterministic model attain very small 
densities of infectives during a significant fraction of the period (i < 10 -5 for 
initial conditions B in the transient regime). It would be interesting to check 
whether the more regular oscillations that have been observed in the standard 
PA for some predator-prey models [7] persist in stochastic simulations of these 
models on RRGs. 
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Fig. 3. For 6 = 1, k = 4, comparison of the solutions of the PA deterministic model 
(blue solid line) with the results of stochastic simulations (black solid line) on a RRG-4 
with N = 10 6 for parameter values in region III, 7 = 0.025, A = 2.5, and two sets of 
initial conditions A and B marked by red dots in the (s, i) plane. Initial conditions: 
a)ssj 0.4889, i « 0.0287, si « 0.0104, sr w 0.2369, ri » 0.0129; b) s R! 0.9240, 
i » 0.0731, si » 0.0558, sr » 0.0024, ri « 0.0005. 




However, the breakdown of the PA SIRS model as the boundary of the oscil- 
latory phase is approached from above and the bistability regime found in region 
III show that there are other effects at play. The standard pair approximation is 
only valid for tree-like structures where each node has exactly the same number 
of contacts and there are no loops, the Bethe lattices. These infinite structures 
cannot be simulated on a computer. On the other hand, classic results of graph 
theory show that a particular realization of a RRG-fc will contain a large number 
of loops, of which the large majority are long (with respect to the average path 
length), so that locally the graph is essentially tree-like. One would expect then 
the PA to perform well on RRGs, provided they are large enough. 

Increasing system size up to N = 10 7 we still find suppression of oscillations 
in region III and significant discrepancies between the transient and steady states 
of the PA SIRS solutions and the results of the simulations in region II close 
to the boundary with region III. A similar problem of oscillation emergence 
and suppression and quantitative differences in Monte Carlo simulations versus 
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mean-field approximation and PA of an evolutionary Rock-Scissors-Paper game 
on different structures was carefully investigated in [8] . For this problem, a more 
accurate multi-site approximation instead of the PA was shown to solve the 
qualitative and quantitative discrepancies with the simulations. The study of 
improved models beyond the PA for SIRS dynamics on RRGs will be the subject 
of future work. 
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